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Abstract 

A new class of methods is introduced for solving the Kohn-Sham equations of density functional 
theory, based on constructing a mapping dynamically between the Kohn-Sham system and an 
auxiliary system. The resulting auxiliary density functional equation is solved implicitly for the 
density response, eliminating the instabilities that arise in conventional techniques for simulations 
of large, metallic or inhomogeneous systems. The auxiliary system is not required to be fermionic, 
and an example bosonic auxiliary density functional is presented which captures the key aspects 
of the fermionic Kohn-Sham behaviour. This bosonic auxiliary scheme is shown to provide good 
performance for a range of bulk materials, and a substantial improvement in the scaling of the 
calculation with system size for a variety of simulation systems. 
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Density functional theory (DFT) [I] has been used successfully to study a vast range of 
chemicals and materials, predicting everything from crystal structures and bulk moduli, to 
XANES and NMR spectra|2j. Its widespread success across the physical sciences has led to 
the creation of entire computational sub-disciplines in holds as diverse as chemistry, physics, 
materials science, engineering and Earth sciences [3], as well as increasing application to 
the biological sciences. Furthermore, DFT simulations have become an essential tool for 
many experimental chemical and materials studies, where they are used to facilitate the 
interpretation of experimental data and guide experimental design. 

In order for a materials simulation to have predictive power, it must take into account the 
quantum mechanical nature of the materials’ constituents, in particular those of the valence 
electrons whose behaviour governs the mechanical, electronic and chemical properties of mat¬ 
ter. In principle the electronic states could be computed by solving the Schrodinger equation 
to find the many-body wavefunction for the electrons, but the computational complexity of 
this approach renders it impractical for realistic systems. DFT provides a solution to this 
problem by showing that only the electronic density is necessary to compute a system’s 
groundstate properties, not the full many-body wavefunction of a system. 

DFT is an exact reformulation of quantum mechanics, founded on the Hohenberg-Kohn 
theorem which states that the groundstate energy of a quantum system is a unique functional 
of the groundstate densityflj[!]. Thus in principle the groundstate energy and density of 
any system may be computed by minimizing the energy with respect to the density, without 
the need to determine the many-body wavefunction. Unfortunately the form of the density 
functional is not known at present. 

Kohn and Sham introduced a practical DFT framework [5J by mapping the ground 
state of the interacting iV-body fermionic system onto that of an auxiliary non-interacting 
fermionic system. This mapping transforms the groundstate solution of the original iV-body 
Schrodinger equation into the N lowest energy solutions of a single-particle Schrodinger 
equation (the Kohn-Sham equation), for which the kinetic, Coulomb and external potential 
functionals have a known, simple form. The Hohenberg-Kohn theorem guarantees that the 
groundstate properties of this auxiliary system will match those of the original many-body 
system, provided they have the same groundstate density; in order for this to be the case, 
an additional functional is included in the auxiliary Kohn-Sham Hamiltonian, referred to as 
the exchange-correlation functional. The exact form of this exchange-correlation functional 
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is not known, but as its contribution to the total energy is small, even relatively crude 
approximations can yield useful results. In this work we extend the mapping of Kohn and 
Sham to a broader class of auxiliary systems, and show how this may be used to improve 
the performance and stability of Kohn-Sham DFT simulations. 

Central to any conventional DFT simulation is the solution of the single-particle Kohn- 
Sham equation. The construction of the charge density is time-consuming and so it is 
common to solve the Kohn-Sham equation for a given input density u in (r): 



where Vh xc [ n] is the Hartree-exchange-correlation potential, V ext is the external potential and 
-0j(r) is the wavefunction of the ith eigenstate (particle) at position r. Clearly in general 
n in (r) is not the same as the density computed from the wavefunctions, n out (r), defined as: 

^out(r) = IVb(r)| 2 , (2) 

3 

where {/-,} are the eigenstate occupation numbers. A separate algorithm is employed to 
evolve u in (r) towards the groundstate density, and it is this algorithm which is the focus of 
this work. Regardless of the method used, once n in (r) has reached the groundstate it must 
satisfy the so-called self-consistency condition 

n in (r) = n out (r) => V Hxc [n in ] = Vh xc [?w]; (3) 

for this reason, these methods are known generally as self-consistent field (SCF) methods. 
Whenever ni n (r) is updated, Equation [l] must be solved to generate the new {^(r)}, and 
hence n ou t(i~)- It is convenient to measure the error in self-consistency at iteration k by the 
density residual R( k \ defined as 

R {k) (r) = (nJS (r) - (r)) • (4) 

When Equations [l] and [ 3 ] are satisfied simultaneously then Ii (k> (r) = 0 and the system has 
reached the ground state. 

The ability of this method to solve the Kohn-Sham equations depends critically on the 
algorithm used to generate nj^(r), the new input density for iteration k, from 1 ' ) (r) and 
i?( fc-1 )( r ). In fact all the common approaches, known as density mixing methods, are based 
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on a linear combination of many previous density residuals, 


n[n 0) = 1} (r) + Z! a mR (m) (r), (5) 

m<k 

where {a m } are real coefficients in the interval [0,1]. 

The popular density mixing methods of PulaypJJ and BroydenJTj work by linearising the 
response of the Kohn-Sham system, assuming that 

n out (r) = M(r) (6) 


where M is a linear operator (though not necessarily Hermitian). In general an approxima¬ 
tion to M is constructed successively by a low-rank update at each iteration, e.g. 


* , ‘> = EEK>M5 1 ("£ ) 

i<kj<k 

where M is a k x k subspace matrix defined as 
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Thus these methods build up an improved approximation to M using the information 
obtained from each iteration. If M is a good representation of the true Kohn-Sham response, 
then the eigendensity of M (as defined by Equation [6| is a good approximation to the self- 
consistent density of the Kohn-Sham system. 

As the size of simulation system is increased, the convergence of simple density mixing 
schemes becomes poorer and these algorithms often diverge for large simulation systems. 
It is well-knownjH] that the root cause of this divergence is the behaviour of the Hartree 
potential Vr, which in reciprocal space may be expressed as 


Vh(G) 


47m in (G) 
O |G| 2 


(9) 


where fi- m is the Fourier transform of n- m , G is a reciprocal lattice vector and is the volume 
of the simulation cell. As the size of simulation cell is increased, the size of the smallest 
G-vector decreases and a small error in n in can lead to a large error Vr [rii n ]; this in turn 
leads to a large error in n out and the iterative scheme will diverge. This instability with 
respect to small changes in n is often referred to as a ‘sloshing instability’. 
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A simple preconditioner was introduced by Kerker |9] in an attempt to correct for this 
ill-behaviour, based on the iterative scheme proposed for jcllium by Manninen et al. ra, and 
this can be incorporated into the density mixing schemes: 

= hhhG) + £ am ( , Jfl 2 ) fl (m, (G) (10) 

m<k VI I ' 0/ 

where Go is a characteristic reciprocal length. The residual components with |G| << Go 
are multiplied by |G| , approximately correcting for the divergence in Vh, whereas the 
components of large wavevectors are unchanged. 

By combining this preconditioner with a density mixing method, reasonable convergence 
can be attained for many moderately sized simulation cells. The characteristic reciprocal 
length Go is a parameter of the scheme, but it is usually found that Go = 1.5A -1 is an 
acceptable approximation for bulk systems [ITT]. 

Preconditioned density mixing methods converge reasonably well for many bulk systems, 
but often struggle with large, metallic or anisotropic simulation systems. There are three 
main reasons for this: firstly, the Kerker preconditioner is isotropic, with a scalar character¬ 
istic reciprocal length Go; secondly the Kerker preconditioner ignores other contributions to 
the response, in particular that due to the exchange-correlation functional; finally the lin¬ 
earisation of the charge density response is only valid near the groundstate, so the methods 
require a good initial input density. 

In principle the exact dielectric response could be calculated from the Kohn-Sham states, 
at least within the linear regime. More practically, the response tensor could be computed 
for the subset containing the smallest G-vectors, a method advocated by Ho et a/.[Hj. Unfor¬ 
tunately this method is unstable when the non-linear response is significant, and furthermore 
the computational cost of the dielectric tensor becomes prohibitive as the simulation size 
increases. This latter issue has been addressed recently by Krotscheck and Liebrecht[l2j, 
who use linear-response perturbation theory to account for the dielectric response. 

The key to making all such schemes tractable is in realising that only a single eigendensity 
is required, and so Equation [6] may be solved efficiently using iterative methods. Such meth¬ 
ods do not require the explicit construction of the response tensor, only the application of 
it to successive trial eigendensities[12]- In this spirit, Raczkowski et a/.|13j proposed replac¬ 
ing the Kerker preconditioner with one based on the Thomas-Fermi-von Weizsacker (TFW) 
functional. In this latter scheme the dielectric response is approximated by a modified TFW 
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system, and the eigendensity computed using an iterative scheme. This eigendensity is then 
used as the input to a conventional Pulay mixing method. 

Whilst this TFW-based scheme avoids the explicit construction of the dielectric tensor, 
it suffers from two major drawbacks: firstly by implementing the method merely as a pre¬ 
conditioner, much of the advantages of the response calculation are lost by the subsequent 
Pulay or Broyden mixing scheme; secondly the iterative preconditioner they proposed is 
not itself stable, and suffers from unphysical run-away solutions. Nevertheless the central 
concept of using an auxiliary functional to model the response is sound. 

In this work a new class of SCF methods is introduced, based on the use of auxiliary den¬ 
sity functionals, which do not require either density-mixing or an explicit dielectric tensor. 
A unified framework for constructing such auxiliary functionals is proposed which avoids the 
separation of the preconditioner and the general diclectic-modelling method and, crucially, 
includes the full non-linear response due to Vh xc - 

In this approach, the Kohn-Sham system is itself modelled using an auxiliary system, 
with a corresponding auxiliary density functional, H 0 [n]. This mapping from the Kohn- 
Sham system to the auxiliary system is exactly analogous to that between the Kohn-Sham 
system and the original many-body interacting system; indeed the existence of an exact 
auxiliary mapping is guaranteed by the Hohenberg-Kohn theorem. 

In keeping with conventional SCF methods, the Kohn-Sham wavefunctions are relaxed 
for a given input density, n in . Once the ground state has been obtained for this input density, 
the Kohn-Sham wavefunctions are used to generate a new density, the output density n out . 
The goal of the SCF method is to find the input charge density which is invariant under 
this process, i.e. to find n i n such that n out = n- m (Equation [3]). As with the Pulay and 
Broyden schemes, the Kohn-Sham response is modelled to provide a good input density 
n m for the next iteration, but rather than linearising the full response and using matrix 
methods, in this approach the Kohn-Sham system is modelled using an auxiliary system, 
with an auxiliary density functional, Ho, which crucially includes the full, non-linear Vh xc - 
This mapping from the Kohn-Sham system to the auxiliary system is exactly analogous to 
that between the Kohn-Sham system and the original many-body interacting system; indeed 
the existence of an exact auxiliary mapping is guaranteed by the Hohenberg-Kohn theorem. 

There are many possible choices for the model functional, each corresponding to a different 
auxiliary system. As an exemplar, we here introduce one of the simplest such choices: a 
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mapping of the fermionic Kohn-Sham system onto an auxiliary system of bosons. This 
bosonic mapping leads to the simple Kohn-Sham-like functional for the auxiliary system 

Ho [n] = (~V 2 + Vhxc [n] + Kxt) • (11) 

A key feature of this auxiliary density functional is that it has the correct behaviour with 
respect to Vhxc and the local V ext , which is required to suppress sloshing instabilities. The 
kinetic and any non-local V ex t contributions are approximations to those of the Kohn-Sham 
system, but are correct in the limit of a simulation system with a single state (the single- 
orbital approximation of Hodgson et al. HI) ; in order for the auxiliary bosonic system to 
model the Kohn-Sham system accurately it is necessary to improve this kinetic-non-local 
approximation dynamically over the course of a simulation. 

Each iteration k of the Kohn-Sham solver (Equation [l| generates a new n'H and n„ut 
pair, which yield information on the non-self-consistent response of the fermionic Kohn-Sham 
system. In order to generate an accurate mapping between the Kohn-Sham and auxiliary 
systems, the auxiliary bosonic system must yield the same density response as the Kohn- 
Sham system when the input density is fixed; since in general this will not be the case, a 
correction operator P is added to the auxiliary Hamiltonian such that at iteration k 



where Ep is the Fermi energy. In order for Equation [12] to be satisfied, 



(13) 


In general P = P [n\, with a corresponding energy correction Ep [n]. This functional is 
unknown, but a Taylor expansion of Ep shows that the first-order variation with n may 
be captured by a fixed operator P. In contrast to density-mixing methods, this linearisa¬ 
tion applies only to the correction energy; in particular the non-linear exchange-correlation 
response is modelled exactly. 

There is considerable freedom in the choice of P, but the success of the Pulay and Broyden 
schemes suggests a low-rank approximation may be productive, for example the Hermitian 
form 


p(fc) 


Qk) (qk 1 
(Qk\Pk) ’ 


(14) 
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(16) 


This form for P ^ is well-defined provided (q k \ pk) ^ 0 is real (which follows naturally 


from the definitions in equations 13, 15 and 16) and non-zero; in the case {q k \ Pk) — 0, 
Equation 13 is satisfied by P = 0. 


Once a suitable auxiliary system has been chosen and the correction operator P de¬ 
termined, the Kohn-Sham system is modelled using the Schrodinger-like equation for the 
square-root of the density 

(17) 


(Ho H + P (k) ) vV = 4‘W 


The (self-consistent) solution to Equation 17 is the ground state density for the auxiliary 


system, n aux , and an estimate of the ground state density of the Kohn-Sham system. This 

density is thus used as the new input density for the Kohn-Sham solver (Equation [I]), i.e. 
(fc+i) 

n- n = n aux . 

As the bosons may all occupy the same state simultaneously, there is no need to compute 
and orthonormalise large numbers of states for the auxiliary system; n (r) is the charge 
density of the entire auxiliary system, with no need for a further summation over bands or 
Brillouin zone sampling points. With an appropriate preconditioner [T5] the computational 


effort to solve Equation [17] scales approximately linearly with system size, and for large 
simulation systems its solution time is insignificant compared to solving Equation [lj despite 
the similarity between the two. 

This auxiliary density functional theory (ADFT) scheme has been implemented in a 
development version of the CASTEP plane-wave pseudopotential computer program[16], 
The charge density is set initially to a superposition of pseudoatomic charge densities, and 
the present scheme compared to the existing Kerker, Pulay and Broyden density mixing 
methods, the latter two methods using a Kerker preconditioner. To illustrate the power of 
the ADFT method, results for a memory-less approximation to P are reported, whereby P^ 
is determined using only information from iteration k; in contrast the Pulay and Broyden 
calculations reported here use information from all previous iterations. 
















Material 


Method 

K A1 Au Si MgO GaN graphite 

Pulay 

6 5 

6 

6 7 8 

6 

ADFT 

5 5 

6 

5 5 6 
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TABLE I: A comparison of the number of SCF iterations required to converge to the ground state 
for a Pulay mixing scheme and the bosonic auxiliary DFT scheme for a variety of simple bulk 
materials. 

The initial tests consisted of applying the bosonic ADFT method to a variety of simple 
materials with different electronic properties, spanning the range from metals to wide band- 
gap insulators and with s-, p- and d-states. Performance was measured by comparing the 
number of iterations required for convergence with that of a conventional Pulay density 
mixing method; in each case the calculation was deemed to have converged when the energy 
was less than 2 peV/atom from the (precomputed) groundstate energy. For all of these tests 
the Pulay and bosonic ADFT methods showed comparable performance (see Table [T]); since 
existing density mixing methods are known to perform well for small, isotropic systems 
such as these the benefit of the ADFT scheme is expected to be for larger, anisotropic 
simulations such as the surface, interface and defect calculations that characterize much of 
nanomaterials, surface science and atomic growth research. 

Magnesium oxide (MgO) is a wide-band-gap insulator with potential applications in 
future CMOS devices. It has been widely studied with ab initio methods, including recent 
studies of electron tunnelling [T7], layer-by-layer film growth [US], and its use as an interlayer 
in stabilising heterostructures|TP]. MgO is also a prototypical polar oxide, and forms an 
ideal testing ground for the application of the ADFT scheme to ionically bonded systems. 

Bulk MgO has a rocksalt structure and conventional density mixing algorithms perform 
well (see Table [[]). In order to simulate thin films the bulk crystal was cleaved along (001) and 
the resulting slabs separated from their periodic image by a vacuum region. The convergence 
behaviour of the ADFT and density mixing methods was then studied as a function of the 
vacuum separation (Figure [I]). 

The strength of the ADFT scheme is apparent as the vacuum gap is increased: the per¬ 
formance of the density mixing schemes degrades considerably, with the number of iterations 
growing approximately linearly with system size, whereas the ADFT performs equally well 
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for all system sizes. The consistent ADFT performance is because the response of the Kohn- 
Sham system is dominated by the Hartree-exchange-correlation and local external potential 
contributions, and these contributions are captured exactly by the auxiliary bosonic system. 

Graphene has a wide variety of potential applications[20], from flexible display tech¬ 
nologies to transistors. It is comprised of a single 2D sheet of carbon hexagons and many 
of its possible applications arise from its unusual bandstructure: graphene is a semi-metal, 
with a Dirac cone at the Fermi level. For this reason the electronic structure of graphene 
and graphene-derivatives are of immense interest, but the computational time required 
by DFT algorithms has led to researchers using quicker tight-binding models for larger 
geometries [21J [22]. 

To simulate graphene a primitive cell was constructed and the convergence behaviour of 
the various methods was studied as the inter-layer separation was increased from 3A (com¬ 
parable to the layer separation in graphite) to 13A by which point each layer is decoupled 
from its periodic images. As the system size increases the iterations required for conventional 
density mixing method grows linearly, whereas the ADFT scheme performs consistently well 
(Figure [2]). 

Gold nanoparticles have a wide variety of medical applications as well as an increasing 
array of applications in electron microscopy and nanomaterials research. Recent DFT sim¬ 
ulations have focused on their high catalytic activity, in particular for oxidation of carbon 
monoxide [ 23 ] • 

To simulate a simple gold nanoparticle a primitive cell was constructed containing a 
gold tetramer, and the inter-particle separation was increased from 7A to 20A at which 
size each particle is decoupled from its periodic images. The convergence behaviour of the 
various methods can be seen in Figure [3] As before, the performance of the ADFT method 
is consistent across the range of cell sizes whereas that of the conventional density mixing 
methods degrades as the vacuum separation is increased, with a linear increase in the number 
of iterations with system size. 

In conclusion, an auxiliary density functional theory framework is introduced that maps 
the fermionic Kohn-Sham system onto an auxiliary system, which is not constrained to be 
fermionic. This mapping is constructed dynamically during an ordinary Kohn-Sham DFT 
calculation, and can be used to accelerate the convergence of SCF calculations considerably. 
The power of this ADFT framework is illustrated using a bosonic auxiliary system and the 
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FIG. 1: (Color online) Convergence of the bosonic ADFT scheme compared to the Pulay, Broyden 
and Kerker (inset) density mixing methods for a 2D slab of MgO as the vacuum gap between 
periodic slab images is varied. 



FIG. 2: (Color online) Performance of the bosonic ADFT scheme for graphene, as a function of the 
spacing between periodic images in the perpendicular direction. The number of iterations taken to 
converge is approximately constant, whereas the iteration count for conventional density mixing 
methods grows. 

performance shown to be invariant with the size of the simulation cell, eliminating the ‘slosh¬ 
ing’ instabilities common in other techniques and removing the linear scaling of iteration 
count with system size. The method has been implemented in a development version of 
the CASTEP computer program, and demonstrates significantly superior performance to 
existing methods for a variety of systems from a wide-band-gap polar oxide thin him, to a 
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FIG. 3: (Color online) Performance of the bosonic ADFT scheme compared to the Pulay, Broyden 
and Kerker (inset) density mixing methods for a gold nanoparticle, as a function of the spacing 
between periodic images. 

2D Dirac semi-metal and an isolated nanocluster. 
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